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ABSTRACT 


A  global-scale  T^-wave  tomography  model  of  the  crust  and  mantle  is  constructed.  The  model  depicts  both  large-scale 
tectonic/dynamic  features,  as  well  as  detailed  upper  mantle  heterogeneities  beneath  the  greater  Middle  East  region. 
Fully  three-dimensional  ray  tracing  is  employed  to  achieve  accurate  travel-time  predictions  of  a  suite  of  P  and  Pn 
arrivals,  necessitating  the  characterization  of  irregular  and  discontinuous  boundaries.  Therefore,  we  explicitly 
represent  undulating  seismic  discontinuities  in  the  crust  and  upper  mantle  within  a  spherical  tessellation  modeling 
framework.  The  tessellation-based  model  architecture  is  hierarchical  in  that  fine  node  sampling  is  achieved  by 
recursively  subdividing  a  base-level  tessellation.  Determining  the  required  node  spacing  to  effectively  model  a  given 
set  of  data  is  problematic,  given  the  uneven  sampling  of  seismic  data  and  the  differing  wavelengths  of  actual  seismic 
heterogeneity.  To  address  this  problem,  we  have  developed  an  inversion  process  called  Progressive  Multi-level 
Tessellation  Inversion  (PMTI)  that  exploits  the  hierarchical  nature  of  the  tessellation-based  design  and  is  well-suited 
for  modeling  a  mixture  of  teleseismic  and  regional  data.  PMTI  allows  the  data  to  determine  the  level  of  model 
complexity  by  progressively  solving  for  shorter  wavelength  structure,  thereby  robustly  imaging  regional  trends  and 
allowing  localized  details  to  emerge  where  resolution  is  sufficient.  Using  the  PMTI  approach,  we  have  developed  a 
global-scale  P-wave  model  that  simultaneously  predicts  teleseismic  P  and  Pn  travel  times  for  events  occurring 
throughout  the  Middle  East.  The  tomographic  image  provides  a  new  glimpse  of  the  complex  upper  mantle  velocity 
anomalies  associated  with  the  convergence  of  the  Arabian  and  Indian  plates  with  Eurasia  along  the  Tethyan  margin. 
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OBJECTIVES 


Travel  time  predictions  based  on  global  tomography  models  can  improve  seismic  event  location  accuracy 
(e.g.,  Antolik  et  al.  2001,  2003).  However,  P-wslyq  models,  derived  from  teleseismic  data,  lack  detailed  structures  in 
the  upper  mantle  and  generally  do  not  accurately  predict  regional  travel  times  as  a  result.  Regional-scale  models  of 
upper  mantle  P-wave  velocity  structure  are  therefore  essential  for  locating  small  seismic  events.  Detailed  upper 
mantle  P-wave  velocities  may  be  interpreted  from  surface  wave  solutions  with  some  success  (e.g.,  Ritzwoller  et  al. 
2003);  however  these  solutions  are  not  specifically  optimized  for  the  prediction  of  Pn  travel  times.  Additionally, 
hybrid  combinations  of  independent  regional  and  teleseismic  models  tend  to  be  inconsistent  since  the  predicted 
teleseismic  (i.e.,  P)  and  regional  (i.e.,  Py)  travel  times  often  conflict  (Yang  et  al.  2004).  At  a  minimum,  baseline 
travel  time  shifts  must  be  determined  through  the  calculation  of  travel  time  correction  surfaces  in  order  to  reduce  the 
D^-order  regional-teleseismic  travel  time  inconsistencies  generated  by  the  hybrid  models. 

The  objective  of  this  project  is  to  further  advance  event  monitoring  capabilities,  through  the  generation  of  a  seamless 
model  of  the  Earth’s  crust  and  mantle  that  is  capable  of  self-consistently  predicting  regional  and  teleseismic  travel 
times  as  accurately  as  possible.  Effectively,  this  requires  simultaneous  modeling  of  detailed  upper  mantle 
heterogeneities  to  adequately  predict  regional  phases  such  as  Pn  and  lower  mantle  structures  to  predict  teleseismic 
travel  times.  It  is  also  highly  desirable  to  develop  a  model  and  techniques  that  reduce  the  need  for  numerous  travel 
time  corrections  (e.g.,  crustal  corrections.  Earth’s  ellipticity,  etc.)  and/or  massive  travel  time  lookup  tables.  In  this 
report,  we  present  our  progress  toward  the  development  of  a  global  3-D  seismic  event  monitoring  model. 

RESEARCH  ACCOMPLISHED 


We  have  established  a  set  of  global  3-D  tomographic  model  design  components  and  procedures  including  a 
statistical  event  location/data  culling  process  and  a  novel  inversion  technique  that  is  specifically  designed  to  exploit 
the  hierarchical  nature  of  our  tessellation-based  model  representation.  We  apply  our  modeling  procedures  to  produce 
a  global-scale  model  of  P-wave  velocity  structure  based  upon  P  and  Pn  travel  times  from  >4,000  Middle  East  events 
and  an  additional  -1,000  globally-distributed  seismic  events.  The  new  tomographic  image  contains  large-scale 
features  seen  in  most  modern  global  tomographic  models,  but  also  contains  previously  un-imaged  details  of  the 
Arabia  plate  beneath  Eurasia.  The  model  is  consistent  with  geology  and  predicts  >800,000  P  and  Pn  travel  times  to 
with  0.53  seconds  standard  deviation  when  statistics-based  data  corrections  are  considered. 

Model  Representation  and  Inversion  Technique 

Our  model  parameterization  is  based  on  triangular  tessellations  of  a  spherical  surface.  Spherical  tessellation 
parameterizations  have  been  employed  for  multiple  geophysical  applications  including  mantle  convection 
simulations  (Baumgardner  and  Frederickson  1985),  magnetic  field  modeling  (Constable  et  al.  1993),  and  the  basis 
for  representation  of  Earth’s  seismic  velocity  structure  (e.g.,  Wang  and  Dahlen  1995;  Chiao  and  Kuo  2001;  Ishii  and 
Dziewonski  2002;  Ballard  et  al.  2009;  Myers  et  al.  2010).  The  primary  purpose  for  employing  spherical  tessellations 
is  to  generate  a  mesh  of  nearly  evenly  spaced  nodes  (or  knots),  therefore  preventing  the  polar  distortion  effects  of  a 
latitude-longitude  grid.  Building  a  spherical  tessellation  mesh  involves  recursive  subdivision  of  triangular  facets  of 
an  initial  polyhedron  (in  our  case,  an  icosahedron)  into  smaller  sub-triangles  whose  vertices  are  projected  onto  a  unit 
sphere  (Figure  1).  Each  sub-division  of  a  triangular  facet  (‘parenf)  produces  4  triangular  facets  {'daughters')  and 
thus  increases  the  number  of  nearly  evenly  spaced  vertices.  We  refer  to  each  tessellation  recursion  as  a  'Level' 
starting  with  the  initial  object  (Level  1;  Figure  1)  up  to  Level  A,  where  N  is  some  maximum  recursion  level.  For 
example,  a  Level  7  recursion  results  in  node  spacing  of  approximately  1  °  when  an  icosahedron  is  used  as  the  initial 
object.  Although  the  tessellation  meshes  represent  simple  sphere-like  surfaces,  they  form  the  basis  of  an  aspherical 
Earth  representation  in  this  study.  See  Ballard  et  al.  (2008,  2009)  and  Simmons  et  al.  (2009)  for  more  information  on 
the  model  design  and  communication. 

The  tomographic  imaging  procedure  developed  for  this  project,  PMTI,  is  designed  with  the  same  underlying 
philosophy  as  the  wavelet  and  multi-grid  approaches  described  in  Chiao  and  Kuo  (2001)  and  Zhou  (1996), 
respectively.  Specifically,  the  PMTI  technique  circumvents  the  need  for  predefining  variable  resolution  scales 
through  the  hierarchical  expression  of  higher  resolution  solutions  as  perturbations  to  lower  resolution  models.  The 
most  notable  difference  of  PMTI  from  the  wavelet  and  multi-grid  approaches  is  that  model  components  are  found 
progressively  through  inversions  of  tomographic  systems  built  upon  successively  finer  parameterizations  defined  by 
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a  hierarchy  of  tessellation  grids.  Therefore,  multiple  parameter  grids  with  different  node  densities  are  not  considered 
simultaneously  and  wavelet  transforms  are  unnecessary.  Each  successive  model  component  is  added  to  the  previous 
solution  to  form  a  total  model  field.  As  with  the  wavelet  and  multi-grid  approaches,  the  data  drive  the  resulting 
variable  resolution  scales  and  finer  detail  structures  emerge  from  the  regional  trends  only  where  it  is  necessary  to 
explain  the  given  set  of  observations.  The  PMTI  process  is  described  at  a  high  level  by  the  following  algorithm: 

x,  =  A;'b 
=b-A^x^ 
for  i  =  2... N 

X,.  =  x*_,  +  A:‘b,._; 
b,=b-A,x, 

end 

Beginning  with  the  initial  travel  time  residuals  (b)  and  sensitivity  kernels  defined  at  the  lowest  resolution  level  (Ay), 
an  initial  slowness  perturbation  model  (xy)  is  determined  through  an  inversion  process.  The  initial  travel  time 
residuals  are  subsequently  reduced  by  the  residuals  predicted  by  the  lowest-resolution  solution  to  produce  a  new 
data  vector,  by.  An  inversion  is  then  performed  on  the  basis  of  sensitivity  kernels  defined  at  the  next  resolution  level 
(A2)  and  the  reduced  data  vector  (by).  The  total  slowness  perturbation  model  (X2)  accumulates  by  summing  the 
interpolated  model  from  the  previous  step  (  and  the  higher  resolution  solution  based  on  the  remaining  travel  time 

signals  (  A“^b^  )•  All  subsequent  inversions  operate  on  data  vectors  reduced  by  the  previous  solution  (b/.y)  to  produce 
slowness  perturbations  relative  to  the  previous  model.  Therefore,  if  a  low-resolution  solution  accounts  for  all  of  the 
travel  time  signal,  the  total  slowness  model  ceases  to  accumulate  additional  slowness  perturbations. 


Figure  1.  Spherical  tessellation  grids  and  parent-daughter  triangle  relationships.  The  initial  object  (a  regular 
icosahedron;  Level  1)  consists  of  20  triangular  faces  and  12  distinct  vertices.  The  triangular  faces 
(parents)  are  recursively  subdivided  into  smaller  triangles  (daughters)  and  the  vertices  are 
normalized  to  the  unit  sphere.  Each  recursion  represents  a  level  in  the  tessellation  hierarchy  and 
parent-daughter  triangle  relationships  are  indexed  for  efficient  communication.  Triangular  regions 
are  shaded  simply  for  visual  orientation  of  the  parent-daughter  regions  and  sub-regions. 

To  evaluate  the  behavior  of  the  PMTI  approach,  we  devised  a  3-D  synthetic  test  bed  to  explore  the  intricacies  of  the 
imaging  procedure.  The  synthetic  model  domain  consists  of  3  triangular  depth  slices  with  nodes  defined  by  a 
tessellation  procedure  up  to  the  Level  5  (Figure  2).  We  chose  to  mimic  complex,  sub-crustal  P  -wave  velocity 
structure  with  two  dominant  characteristic  wavelengths.  Specifically,  the  synthetic  model  consists  of  a  long- 
wavelength  regional  velocity  trend  with  short- wavelength  details  embedded  within  the  center  of  model  domain 
(Figure  2,  left).  The  chosen  ray  path  configuration  consisted  of  clusters  of  regional  phases  (sub-lateral  travel 
direction)  and  teleseismic  paths  (near  vertical  travel  direction).  The  ray  path  density  is  highly  heterogeneous  and 
random  noise  was  added  to  imitate  a  realistic  shallow  mantle  imaging  problem. 

Inversions  were  performed  directly  (simple  least-squares  at  the  highest  tessellation  level)  and  with  the  PMTI 
procedure  for  comparison.  A  representative  solution  from  a  single  realization  (paths  shown  in  Figure  2  and  20% 
RMS  noise  added)  is  shown  in  Figure  3.  For  the  simple  (or  direct)  inversion,  we  incorporated  damping  and 
smoothing  operators  independently  and  optimized  the  regularization  weights  to  find  the  solutions  closest  to  the 
known  solution  (based  on  RMS  model  misfit).  With  damping  alone,  detailed  structures  emerge  where  there  are  a 
large  number  of  crossing  paths  (Figure  3;  SIMPLE  damping).  However,  instead  of  recovering  the  long-wavelength 
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regional  trend,  isolated  velocity  structures  are  generated  without  a  smooth  transition  through  unresolved  regions.  In 
regions  with  only  near-vertical  paths  (teleseismic  phases),  circular  velocity  structures  are  generated  since  there  is  no 
sensitivity  to  surrounding  nodes.  In  regions  with  only  bundles  of  sub-horizontal  paths  (regional  phases),  linear 
velocity  anomalies  are  evident.  Such  isolated  velocity  blobs  and  streaks  are  undesirable  since  they  may  result  in 
geological  misinterpretations  and  also  pose  serious  problems  when  computing  travel  times  with  3-D  ray  paths. 
Alternatively,  simple  inversion  with  a  first-order  smoothing  operator  (with  an  optimum  weight)  recovers  the  long- 
wavelength  regional  trend,  but  fails  to  recover  fine  details  where  information  is  sufficient  (Figure  7;  SIMPLE 
smoothing).  This  smooth  solution  is  also  undesirable  since  we  would  prefer  to  take  full  advantage  of  the 
resolvability  of  dense  data  packets. 


B  B2  e.fl 

Vp  (kiTi/S) 


Figure  2.  Synthetic  tomography  scenario  imitating  a  combination  of  regional  and  teleseismic  phases  and 

complex  shallow  upper  mantle  velocity  structure  parameterized  by  a  triangular  tessellation  with  5 
hierarchical  levels,  (left)  Synthesized  velocity  structure  with  a  regional  trend  (from  left  to  right)  and 
embedded,  small-scale  velocity  heterogeneities,  (center)  Map  view  of  synthetic  ray  path  coverage 
consisting  of  several  clusters  of  sub-horizontal  regional  phase  ray  paths  and  3  clusters  of  near¬ 
vertical  teleseismic  ray  paths  (highlighted  with  green  circles),  (right)  3-D  view  of  the  synthesized  ray 
path  coverage. 

Using  the  PMTI  approach,  the  long-wavelength  velocity  structure  is  established  with  the  first  inversion  based  on 
nodes  defined  by  the  first  level  tessellation  (Figure  3;  PMTI  Level  1).  Note  that  only  simple  damping  is  used  as 
regularization.  The  Levels  2-4  PMTI  solutions  adjust  the  model  only  slightly,  relative  to  the  Level  1  solution.  The 
minor  model  variation  provided  by  these  intermediate  steps  is  due  to  the  synthetic  velocity  model  design  with  only 
two  dominant  characteristic  wavelengths.  In  the  final  PMTI  stage  (Level  5  inversion),  detailed  structures  emerge 
from  the  underlying  long- wavelength  velocity  trend  (Figure  3;  PMTI  Level  5).  The  detailed  anomalies  are  most 
refined  where  data  coverage  is  dense  and  ray  path  travel  directions  are  diverse.  It  should  be  noted  that  the  slight 
changes  through  the  first  3  progressions  demonstrates  the  potential  pitfalls  of  an  alternative  imaging  procedure 
involving  dynamic  mesh  refinement  (e.g.,  Sambridge  and  Faletic  2003).  A  dynamic  mesh  refinement  procedure 
would  likely  cease  to  insert  new  model  nodes,  thereby  never  recovering  the  fine  (yet  resolvable)  details.  Our  tests 
also  reveal  that  the  PMTI  solution  is  less  dependent  on  the  choice  of  damping  weight  than  a  simple  direct  inversion. 
While  it  is  true  that  a  similar  model  could  be  recovered  through  a  simple  inversion  with  combined  damping  and 
smoothing,  regularization  weights  would  need  to  be  fine-tuned  in  order  to  best  match  the  known  solution.  In  a  real 
scenario,  however,  the  ‘true’  solution  is  never  known. 

The  PMTI  approach  has  important  advantages  over  other  inversion  techniques,  especially  when  considering  extreme 
spatial  resolvability  variations.  If  the  travel  time  observations  cannot  be  explained  by  a  long- wavelength  velocity 
trend  in  a  particular  region,  detailed  anomalies  emerge  from  the  background  velocity  model.  Therefore,  the 
‘resolution  level’  is  predicated  by  the  data  without  any  user  intervention.  Dynamic  mesh  refinement  (or  adaptive 
meshing)  is  one  such  method  of  intervention  whereby  the  addition  or  subtraction  of  model  parameters  occurs  based 
on  some  user-defined  criteria,  thus  generating  an  irregular  grid.  Aside  from  the  additional  complexity  of  establishing 
and  communicating  with  a  model  based  on  an  irregular  grid,  our  synthetic  example  demonstrates  other  potential 
pitfalls  of  a  dynamic/adaptive  mesh  refinement  approach  (Figure  3).  For  example,  if  we  chose  to  refine  the  mesh 
from  one  inversion  step  to  the  next  based  on  model  movement,  the  imaging  process  would  cease  after  the  inversion 
based  on  the  nodes  defined  by  the  Level  2  triangular  tessellation  (see  Figure  3).  This  cessation  could  be  prevented 
through  knowledge  of  the  data  density  (e.g.,  hit  count  criteria).  However,  hit  count  does  not  necessarily  reflect 
resolvability  since  ray  paths  are  commonly  bundled  providing  limited  additional  constraints  on  lateral  velocity 
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variations.  For  ultra-high  resolution  modeling  with  a  global  basis  (<0.5°  node  spacing),  adaptive  meshing  may  be 
required  to  limit  the  size  of  the  tomographic  systems  of  equations.  For  our  initial  global  model  design  however,  we 
find  that  irregular  meshing  is  unnecessary. 


-500  0  500  -500  0  500  -SOO  0  500  -500  0  500 


X(km)  X(km)  X(km)  X(lcm) 

Figure  3.  Representative  tomography  solutions  for  the  synthetic  tomography  scenario  shown  in  Figure  2.  (top 
row)  Synthetic  model  and  3  possible  solutions  using  different  approaches  and  regularization.  With 
simple  inversion  and  damping  (^SIMPLE  damping’),  detailed  structures  emerge  where  there  are 
sufficient  crossing  paths  and  isolated  velocity  structures  are  generated  where  there  are  gaps  in  ray 
path  coverage.  Simple  inversion  with  a  first-order  smoothing  operator  (^SIMPLE  smoothing’) 
recovers  the  long-wavelength  regional  trend,  but  fails  to  recover  fine  details.  Using  the  PMTI 
approach  (TMTI  Level  5’),  the  regional  trend  is  recovered  and  details  emerge  where  there  is 
sufficient  ray  path  coverage  (in  the  northern  portion  of  the  model  space),  (bottom  row)  The 
accumulation  of  velocity  structure  through  the  PMTI  process. 

Travel  Time  Data  and  Event  Locations 

We  compiled  a  Middle  East-centric  set  of  teleseismic  P  and  regional  travel  time  measurements  from  a  subset  of 
our  local  database  at  Lawrence  Livermore  National  Laboratory  (Ruppert  et  al.  2005)  and  the  publicly  available 
Engdahl-van  der  Hilst-Buland  (EHB)  catalog  (Engdahl  et  al.  1998).  The  data  selection  criteria  were:  1)  all  available 
P  and  Pn  travel  time  measurements  from  seismic  events  occurring  within  the  greater  Middle  East  region,  and  2) 
teleseismic  P  measurements  from  globally  distributed  events  with  the  largest  number  of  recordings  within  5°  lateral 
bins  and  7  depth  bins  from  the  surface  to  700  km  depth.  In  total,  we  consider  5,401  events  (~4,000  within  the 
Middle  East)  recorded  at  approximately  4,500  seismic  stations  around  the  globe.  The  data  consist  of  -800,000  P 
arrivals  from  globally  distributed  events  and  stations,  and  -43,000  Pn  arrivals  to  stations  within  the  Middle  East 
(Figure  4).  Although  the  compiled  dataset  is  global  in  extent,  more  than  60%  of  the  data  are  from  events  occurring 
within  the  Middle  East  and/or  arrivals  to  one  of  the  600+  seismic  stations  within  the  region.  Ray  path  hitcount 
calculations  reveal  that  most  of  the  shallow  upper  mantle  in  the  Middle  East  region  is  sampled  by  more  than  1,000  P 
and  100  Pn  paths  providing  excellent  constraints  on  detailed  velocity  structure  (Figure  4). 

We  adapted  a  multi-event  location  algorithm,  initially  designed  to  simultaneously  locate  event  clusters,  to  handle  a 
global  distribution  of  seismic  events  (Myers  et  al.  2007,  2009).  The  algorithm,  called  BayesLoc,  is  a  multi-event 
location  algorithm  that  operates  within  a  hierarchical  Bayesian  statistical  framework.  The  process  entails  generating 
numerous  realizations  of  event  hypocenters,  origin  times,  phase  labels,  and  travel  time  curves  using  a  Markov  Chain 
Monte  Carlo  (MCMC)  sampling  scheme.  Using  prior  statistical  models,  the  process  allows  for  complete  statistical 
characterization  of  the  multi-event  problem  and  thus  develops  a  joint  posterior  distribution  of  all  elements  involved. 
The  BayesLoc  algorithm  simultaneously  determines  event  location  probability  regions,  identifies  likely  outlier 
arrival  time  data  (inconsistent  with  the  total  population),  re-labels  misidentified  seismic  phases,  and  adjusts  travel 
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times  based  on  statistically  likelihood.  The  net  result  of  the  procedure  is  a  set  of  event  locations  and  travel  times 
such  that  the  spread  of  the  travel  time  residuals  relative  to  the  1-D  AK135  (Kennett  et  al.  1995)  velocity  model  is 
minimal  (Figure  5). 


In  terms  of  data  variance,  the  apparent  travel  time  residual  signal  of  the  raw  data  (based  on  the  original  event 
locations)  is  ~12  times  larger  than  the  signal  based  on  the  BayesLoc  event  locations.  This  suggests  that  a  tremendous 
amount  of  the  apparent  signal  associated  with  the  raw  event  locations  is  not  due  to  3-D  velocity  structure,  but  rather 
hypocenter  mis-location,  mis-identified  seismic  phases  and  arrival  time  pick  errors.  The  total  signal  reduction 
signifies  the  improved  data  consistency  that  BayesLoc  provides,  but  another  important  aspect  is  the  shift  of  the  mean 
from  a  positive  0.44  seconds  to  a  negative  0.53  seconds.  This  mean  shift  and  the  residual  trend  with  distance 
(see  Figure  5)  have  tremendous  implications  on  the  absolute  velocity  structure  needed  to  explain  the  data. 
Hypocenter  locations  based  on  single-event  location  procedures  often  do  not  recover  such  trends;  particularly  when 
assuming  a  location  such  that  the  residual  travel  times  are  zero-mean. 


Stations 


Figure  4.  (left  column)  Selected  seismic  events  (5,401)  and  stations  (4,500)  around  the  globe,  (right  column) 
Pn  and  teleseismic  P  hitcount  in  the  outlined  region  at  a  depth  of  150  km. 
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Figure  5.  Travel  time  residuals  for  the  compiled  P  and  Pn  data  set.  Residual  travel  times  are  plotted  as  a 

function  of  great  circle  distance  in  terms  of  point  density  on  a  normalized  log  scale,  (a)  Travel  time 
residuals  for  the  initial  EHB  bulletin  and  LLNL  data  set  with  respect  to  (w.r.t.)  the  AK135  1-D 
velocity  model,  (b)  Travel  time  residuals  w.r.t.  AK135  after  BayesLoc  processing,  (c)  Residuals 
based  on  BayesLoc  processed  data  relative  to  the  hybrid  regional-teleseismic  velocity  model 
described  in  the  text  (RSTT+GyPSuM).  (d)  Travel  time  residuals  w.r.t.  to  the  global  3-D  P  wave 
model  constructed  in  this  study  (LLNL-G3Dvl.O). 
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3-D  Ray  Tracing  and  Sensitivity  Definitions 

We  have  adapted  a  3-D  ray  tracing  procedure  specifically  tailored  to  communicate  with  the  aforementioned 
tessellation-based  model  design.  The  general  approach  we  have  taken  is  based  on  the  techniques  originally 
conceived  by  Zhao  et  al.  (1992)  to  model  direct  arrivals  within  a  subduction  zone  and  extended  to  teleseismic  phases 
by  Zhao  and  Lei  (2004).  The  procedure  combines  pseudo-bending  in  continuous  media  (Um  and  Thurber  1987) 
while  simultaneously  satisfying  Snell’s  Law  at  discontinuous  interfaces  through  an  iterative  process.  One  major 
shortcoming  of  this  general  approach  is  the  inability  to  find  global  minimum  travel  times  and  paths  for  complex 
regional  phases  such  as  We  have  therefore  adapted  the  procedure  to  find  such  global  minima  by  considering 
multiple  starting  ray  path  configurations.  A  similar  adaptation  has  recently  been  developed  by  Ballard  et  al.  (2009) 
where  the  computational  process  was  designed  within  a  distributed  computing  environment. 

In  our  specific  approach,  we  initially  create  several  crude  3-D  ray  path  estimates  based  on  an  imposed  set  of  rules 
for  a  given  source-receiver  configuration.  Each  of  the  initial  paths  is  perturbed  with  a  limited  set  of  pseudo-bending 
iterations  and  discontinuity  piercing  point  adjustments  as  described  in  Zhao  et  al.  (1992).  Rather  than  performing  ray 
path  adjustment  iterations  until  convergence,  we  cull  out  paths  that  provide  travel  times  within  some  tolerance  of  the 
minimum  time  achieved  by  the  full  set  of  perturbed  paths  (we  chose  a  1  second  tolerance).  Ray  path  adjustment 
iterations  are  then  performed  on  the  remaining  paths  until  the  minimum  travel  time  is  attained.  The  paths  that 
generate  the  minimum  time  within  a  secondary  time  tolerance  (we  chose  0.2  seconds)  are  accepted  and  used  in  the 
development  of  model  space  sensitivity  kernels.  For  example,  near  the  upper  mantle  transition  zone  crossover 
distance  range  (-18°  distance  for  a  surface  event),  ray  paths  traveling  through  the  shallow  upper  mantle  and  below 
transition  zone  discontinuities  are  accepted  if  they  satisfy  the  minimum  time  tolerance  criterion.  By  accepting  all 
paths  that  produce  near-minimum  travel  times,  we  directly  address  the  multi-pathing  problem.  Additionally,  model 
sensitivity  may  be  distributed  over  a  wide  depth  range  as  opposed  to  two  (or  more)  distinct  geological  units  in  a  pure 
multi-pathing  scenario  (see  the  example  in  Figure  6).  Thus,  the  issue  of  interdependence  of  velocity  structure  and 
ray  path  is  mitigated. 


Distance  (degrees) 

Figure  6.  Example  3-D  ray  paths  calculated  through  the  global  velocity  model  constructed  in  this  study. 

Sensitivity  is  spread  across  broad  depth  zones  and/or  multiple  model  units  to  mitigate  the  problem 
of  interdependence  of  paths  and  velocity  structure  as  well  as  multi-pathing. 

Inversion  Procedures  and  Results 

We  designed  a  starting  model  by  combining  the  Regional  Seismic  Travel  Time  (RSTT)  model  (Myers  et  al.  2010) 
and  the  global  joint  seismic-geodynamic  image  called  GyPSuM  (Simmons  et  al.  2010).  The  crust  and  upper  mantle 
velocity  model  is  represented  by  33  layers  (including  7  crustal  units)  with  a  lateral  spacing  of-l°  (Level  7 
tessellation)  and  the  lower  mantle  is  represented  by  25  layers  with  a  lateral  spacing  of  -2°  (Level  6  tessellation).  In 
total,  the  crust  and  mantle  T^-wave  velocity  structure  is  defined  at  -1.6  million  points  (1.35  million  points  in  the  crust 
and  upper  mantle;  250,000  points  in  the  lower  mantle).  Sensitivities  to  each  of  the  model  nodes  defined  in  the 
starting  solution  (RSTT+GyPSuM)  were  determined  based  on  the  calculated  3-D  ray  paths  and  model  space 
sensitivity  definitions  described  in  the  previous  section.  For  the  inverse  problem,  we  chose  to  combine  certain  layers 
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to  limit  the  number  of  free  parameters  and  account  for  limited  resolvability  in  particular  depth  zones.  Specifically, 
sensitivities  to  all  nodes  representing  the  crust  (as  defined  in  the  full  hybrid  model)  were  summed  along  the 
geocentric  vertices  to  define  a  single  set  of  sensitivity  kernels  associated  with  the  crust.  Therefore,  tomographic 
inversion  results  in  slowness  perturbation  updates  for  the  entire  crustal  stack  at  each  lateral  position  instead  of 
determining  model  updates  within  each  individual  crustal  layer.  The  upper  mantle  sensitivities  were  condensed  in  a 
similar  fashion  creating  9  total  inversion  layers  in  the  crust  and  upper  mantle  with  an  average  depth  spacing  of 
~80  km.  The  lower  mantle  sensitivities  were  only  slightly  condensed  providing  an  average  depth  spacing  of  just 
over  100  km.  In  total,  the  node  sensitivities  defined  for  the  starting  model  were  condensed  from  55  model  layers  and 
-1.6  million  nodes  to  31  layers  and  -600,000  free  parameters  to  be  determined  in  the  tomographic  inversion 
process. 

Global-scale  velocity  model  updates  were  determined  by  carrying  out  the  PMTI  procedure  where  successive 
inversions  were  performed  using  the  LSQR  algorithm  (Paige  and  Saunders  1984).  We  initiated  the  progression  by 
solving  for  model  updates  at  nodes  defined  by  tessellation  Level  4  vertices  (-8°  lateral  node  spacing)  followed  by 
Levels  5,6  and  7  (4°,  2°  and  1°  node  spacing,  respectively).  We  chose  not  to  begin  the  progressive  inversion  below 
Level  4  to  limit  structural  leakage  across  drastically  different  tectonic  environments.  In  addition,  we  chose  not  to 
adjust  the  depth  resolution  during  the  progression.  A  global  damping  operator  was  employed  as  regularization  and 
the  damping  weight  was  chosen  based  on  the  trade-off  of  model  norm  and  data  misfit  (L-curve  analysis). 

Owing  to  the  substantial  non-linearity  of  ray  paths  and  velocity  structure  for  regional  phases,  Pn  ray  paths  and 
sensitivity  kernels  were  re-computed  after  one  complete  PMTI  cycle  and  the  inversion  procedure  was  performed 
again.  After  the  2^^  inversion  cycle,  Pn  ray  paths  were  computed  once  again  to  test  for  convergence.  It  was 
concluded  that  most  minimum-time  ray  paths  did  not  change  appreciably;  and  paths  that  noticeably  changed,  did  not 
generate  substantially  different  travel  time  predictions  relative  to  the  previous  stage.  It  is  our  conclusion  that  the 
sensitivity  definitions,  which  include  multi -pathing  and  broadened  kernels,  mitigate  the  non-linear  path- velocity 
interdependence  issues  to  a  large  degree.  However,  we  are  cautious  about  generalizing  this  observation  to  every 
scenario  and  region. 

The  final  P  -wave  model,  LLNL-G3DvL0,  predicts  the  BayesLoc  posteriori  data  to  within  0.53  seconds  standard 
deviation  (Figure  5).  This  equates  to  -72%  variance  reduction  relative  to  BayesLoc  posteriori  travel  times  predicted 
with  AK135  (1.01  second  standard  deviation).  The  initial  bulletin  data  {yNithout  BayesLoc  statistical  operations)  are 
predicted  to  within  3.45  seconds  standard  deviation  on  the  basis  of  the  1-D  AK135  Earth  model.  The  combination  of 
BayesLoc  processing  and  the  LLNL-G3DvL0  velocity  model  yields  data  fits  on  the  order  of  97%  variance  reduction 
relative  to  the  initial  data  and  AK135  (see  Figure  5).  Thus,  the  travel  times  for  the  -800,000  P  and  -43,000  Pn 
arrivals  are  predicted  to  extremely  high  degrees  given  the  statistical  corrections  provided  by  BayesLoc  and  the 
3-D  velocity  heterogeneities  in  the  new  global-scale  tomographic  model  shown  in  Figures  7  and  8. 
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Figure  7.  Selected  360°  cross  sections  through  the  LLNL-G3Dvl.O  velocity  model.  Most  of  the  major  features 
highlighted  are  also  observed  in  the  hybrid  starting  model  based  on  the  joint  seismic-geodynamic 
image  called  GyPSuM  (Simmons  et  al.  2010).  However,  a  number  of  additional  features  are 
observed  in  the  LLNL-G3Dvl.O  model  including  deep  slabs  originating  from  the  Japan  subduction 
zone  and  detailed  upper  mantle  structures  beneath  the  Middle  East  region  (see  Figure  8). 
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Figure  8.  (left  column)  Close-up  view  of  the  P  -wave  velocity  at  150  km  depth  from  the  LLNL-G3Dvl.O  model. 
The  final  model  differs  dramatically  from  the  starting  solution  in  the  shallow  upper  mantle 
beneath  the  Middle  East  region.  Most  notably,  several  high-velocity  structures  emerge  beneath  the 
region  and  correlate  well  with  the  location  of  deep  earthquakes  used  in  this  study,  (right  column) 
Cross  sections  along  the  A-A’,  B-B’  and  C-C’  transects.  Velocity  perturbations  are  relative  to  the 
global  mean  velocity  anomalies  in  LLNL-G3Dvl.O  at  each  depth.  High-velocity  structures 
extending  from  Arabia  to  beneath  Eurasia  represent  underthrusted  continental  lithosphere  along 
the  western  plate  boundary  and  subducted  oceanic  lithosphere  beneath  the  Makran  region  along 
the  southern  plate  boundary. 


CONCLUSIONS  AND  RECOMMENDATIONS 


Our  initial  monitoring  model,  named  LLNL-G3Dvl.O,  is  a  global-scale  tomographic  model  of  P  -wave  velocity  that 
predicts  >800,000  P  and  >46,000  Pn  travel  times  (mostly  from  Middle  East  events)  to  within  0.53  seconds  standard 
deviation.  The  model  is  built  with  the  tools  and  processes  described  in  this  report  including:  BayesLoc  data 
processing,  3-D  ray  tracing  with  multi-pathing,  and  PMTI  imaging.  The  current  model  consists  of  large-scale  mantle 
structures  such  as  spreading  centers,  cratons,  superplumes  and  lower  mantle  slabs  (generally  observed  in  all  modern 
global  tomography  models).  In  addition,  the  image  depicts  remarkably  detailed  structures  in  the  upper  mantle  below 
the  Middle  East,  owing  to  the  dense  data  coverage  in  the  region  and  data  consistency  yielded  by  the  BayesLoc 
processing.  In  addition  to  predicting  travel  times,  our  model  provides  a  new  and  important  perspective  of  the  upper 
mantle  structures  associated  with  convergence  of  the  Arabia  and  India  plates  with  Eurasia.  In  particular,  we  image 
the  underthrusted  continental  Arabian  lithosphere  beneath  Iran  along  the  western  Arabia-Eurasia  plate  boundary.  In 
addition,  the  transition  from  continent-continent  collision  to  subduction  beneath  the  Makran  region  along  the 
southern  Arabia-Eurasia  boundary  is  evidenced  by  the  imaged  structures  in  the  LLNL-G3DvL0  model.  The  next 
phase  in  this  modeling  effort  will  be  to  perform  the  same  procedures  described  in  the  current  study  with  a  massive 
set  of  teleseismic  P  travel  time  data  that  are  more  globally  distributed  than  the  data  considered  herein.  In 
conjunction,  P^  travel  time  data  from  multiple  regions  will  be  considered  simultaneously  for  the  purpose  of  self- 
consistent  travel  time  prediction  at  both  regional  and  teleseismic  distances. 
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